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Abstract 

A practical computational scheme based on time-dependent density functional theory (TDDFT) 
and ultrasoft pseudopotential (USPP) is developed to study electron dynamics in real time. A 
modified Crank-Nicolson time-stepping algorithm is adopted, under planewave basis. The scheme 
is validated by calculating the optical absorption spectra for sodium dimer and benzene molecule. 
As an application of this USPP-TDDFT formalism, we compute the time evolution of a test electron 
packet at the Fermi energy of the left metallic lead crossing a benzene-(l,4)-dithiolate junction. A 
transmission probability of 5-7%, corresponding to a conductance of 4.0-5.6 /iS, is obtained. These 
results are consistent with complex band structure estimates, and Green's function calculation 
results at small bias voltages. 

PACS numbers: 71.15.-m, 73.63.-b, 78.67.-ii 
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I. INTRODUCTION 



The development of molecular scale electronic devices has attracted a great deal of interest 
in the past decade, although major experimental and theoretical challenges still exist.— 
To date precise experimental control of molecular conformation is lacking, resulting in large 
uncertainties in the measured conductance. On the theory side, while the Green's function 
(GF) method has achieved many successes in describing electron transport at the mesc®'^ and 
molecula r^i^i^°i^^i^^ scales, issues such as dynamical electron correlation and large electron- 
phonon coupling effectsi^^i^ are far from fully resolved. It is therefore desirable to exploit 
alternative approaches for comparison with the mainstream GF calculations. In this paper, 
we describe a first step towards this goal by computing how an electron propagates through 
a molecular junction in real time, based on the time-dependent density functional theoryi^ 
(TDDFT). 

Density functional theory (DFT)^^ with the Kohn-Sham reference kinetic energy func- 
tional of a fictitious non-interacting electron system^^ is a leading method for treating many 
electrons in solids and molecules.—. While initially formulated to describe only the elec- 
tronic ground statei^*i^, it has been rigorously extended by Runge and Gross^^ to treat 
time-dependent, driven systems (excited states). TDDFT is therefore a natural theoreti- 
cal platform for studying electron conduction at the nanoscale. There are two flavors in 
which TDDFT is implemented. One is direct numerical integratio n-'^^i^^i^^i^^i^^i^'^ of the time- 
dependent Kohn-Sham (TDKS) equations. The other is a Gedanken experiment of the 
former with an added assumption of infinitesimal time-dependent perturbation, so a linear 
response function may be first derived in closed for m^^i^^i^^ , which is then evaluated numer- 
ically. These two implementations should give exactly the same result when the external 
perturbation field is infinitesimal. The latter implementation can be computationally more 
efficient once the linear-response function has been analytically derived, while the former 
can treat non-infinitesimal perturbations and arbitrary initial states. 

A key step of the TDDFT dynamics is updating of the Kohn-Sham effective potential 
by the present excited-state charge density p(x, t), VKs(t) = Vks[p(x, t), ...]. This is what 
sets TDDFT apart from the ground-state DFT estimate of excitation energies, even when 
TDDFT is applied in its crudest, so-called adiabatic approximation,^^ whereby the same 
exchange-correlation density functional form as the ground-state DFT calculation is used 
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(for example, the so-called TDLDA approximation uses exactly the same Ceperley-Alder- 
Perdew-Zunger functiona P^i^^ as the ground-state LDA calculation.) This difference in ex- 
citation energies comes about because in a ground-state DFT calculation, a virtual orbital 
such as LUMO (lowest unoccupied molecular orbital) experiences an effective potential due 
to electrons occupying the lowest A^ orbitals; whereas in a TDDFT calculation, if one 
electron is excited to a LUMO-like orbital, it sees A^ — 1 electrons occupying the lowest 
A^ — 1 orbitals, plus its own charge density. Also, the excitation energy is defined by the 
collective reaction of this coupled dynamical system to time-dependent perturbation (pole 
in the response function)^, rather than simple algebraic differences between present vir- 
tual and occupied orbital energies. For rather involved reasons beyond what is discussed 
here, TDDFT under the adiabatic approximation gives significantly improved excitation 
spectr a^^i^^ , although there are still much to be desired. Further systematic improvements 
to TDDFT such as current density functional^ and self-interaction correction-^ have already 
made great strides. 

Presently, most electronic conductance calculations based on the Landauer transmission 
formalism!^^^ have assumed a static molecular geometry. In the Landauer picture, dis- 
sipation of the conducting electron energy is assumed to take place in the metallic leads 
(electron reservoirs), not in the narrow molecular junction (channel) itself.— Inelastic scat- 
tering, however, does occur in the molecular junctions themselves, the effects appearing as 
peaks or dips in the measured inelastic electron tunneling spectra (lETS)^ at molecular 
vibrational eigen-frequencies. Since heating is always an important concern for high-density 
electronics, and because molecular junctions tend to be mechanically more fragile compared 
to larger, semiconductor-based devices, the issue of electron-phonon coupling warrants de- 
tailed calculationa^^*^ (here we use the word phonon to denote general vibrations when 
there is no translational symmetry). In the case of long vr-conjugated polymer chain junc- 
tions, strong electron-phonon coupling may even lead to new elementary excitations and 
spin or charge carriers, called soliton/polaroni^»^^'^2»^24ii^ where the electronic excitation is 
so entangled with phonon excitation that separation is no longer possible. 

In view of the above background, there is a need for efficient TDDFT implementations 
that can treat complex electron-electron and electron-phonon interactions in the time do- 
main. Linear-response type analytic derivations can become very cumbersome, and for some 
problems^ may be entirely infeasible. A direct time-stepping metho d^'^'^' - 'i^'^i^^i^^'^i^^ analogous 
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to molecular dynamics for electrons as well as ions may be more flexible and intuitive in 
treating some of these highly complex and coupled problems, if the computational costs can 
be managed. Such a direct time-stepping code also can be used to double-check the correct- 
ness of analytic approaches such as the non-equilibrium Green's function (NEGF) method 
and electron-phonon scattering calculation o^^i^^ , most of which explicitly or implicitly use 
the same set of TDDFT approximations (most often an adiabatic approximation such as 
TDLDA). 

Two issues are of utmost importance when it comes to computational cost: choice 
of basis and pseudopotential. For ground-state DFT calculations that involve a signifi- 
cant number of metal atoms (e.g. surface catalysis), the method that tends to achieve 
the best cost-performance compromise is the ultrasoft pseudopotential (USPP)^-2ii2iii with 
planewave basis, and an independent and theoretically more rigorous formulation, the projec- 
tor augmented-wave (PAW)^ method. Compared to the more traditional norm- conserving 
pseudopotential approaches, USPP/PAW achieve dramatic cost savings for first-row p- and 
(i-elements, with minimal loss of accuracy. USPP/PAW are the workhorses in popular codes 
such as VASP— and DACAPO ^^i^^i^^ . We note that similar to surface catalysis problems, 
metal- molecule interaction at contact is the key for electron conduction across molecular 
junctions. Therefore it seems reasonable to explore how TDDFT, specifically TDKS under 
the adiabatic approximation, performs in the USPP/PAW framework, which may achieve 
similar cost-performance benefits. This is the main distinction between our approach and 
the software package Octopu o^^'^^ , a ground-breaking TDDFT program with direct time 
stepping, but which uses norm-conserving Troullier-Martins (TM) pseudopotential^, and 
real-space grids. We will address the theoretical formulation of TD-USPP (TD-PAW) in 
sec. II, and the numerical implementation of TD-USPP in the direct time-stepping flavor 
in sec. III. 

To validate that the direct time-integration USPP-TDDFT algorithm indeed works, we 
calculate the optical absorption spectra of sodium dimer and benzene molecule in sec. IV and 
compare them with experimental results and other TDLDA calculations. As an application, 
we perform a computer experiment in sec. V which is a verbatim implementation of the 
original Landauer pictur o^^i'^^ . An electron wave pack comes from the left metallic lead (ID 
Au chain) with an energy that is exactly the Fermi energy of the metal (the Fermi electron), 
and undergoes scattering by the molecular junction (benzene- (l,4)-dithiolate, or BDT). The 
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probability of electron transmission is carefully analyzed in density vs. x, t plots. The 
point of this exercise is to check the stability and accuracy of the time integrator, rather 
than to obtain new results about the Au-BDT-Au junction conductance. We check the 
transmission probability thus obtained with simple estimate from complex band structure 
calculationsSiiSS,^ and Green's function calculations at small bias voltages. Both seem to be 
consistent with our calculations. Lastly, we give a brief summary in sec. VI. 

II. TDDFT FORMALISM WITH ULTRASOFT PSEUDOPOTENTIAL 

The key idea of USPP/PAW— i^^*^^*^ is a mapping of the true valence electron wave- 
function ip{yi) to a pseudowavefunction ipi^)'- 'ip ^ i', like in any pseudopotential scheme. 
However, by discarding the requirement that ipi^) must be norm-conserved {{iplip) = 1) 
while matching outside the pseudopotential cutoff, a greater smoothness of ^/'(x) in the 
core region can be achieved; and therefore less planewaves are required to represent '?/'(x). In 
order for the physics to still work, one must define augmentation charges in the core region, 
and solve a generalized eigenvalue problem 

H\^ljn)=eJ\^n), (1) 

instead of the traditional eigenvalue problem, where S" is a Hermitian and positive definite 
operator. S specifies the fundamental measure of the linear Hilbert space of pseudowave- 
functions. Physically meaningful inner product between two pseudowavefunctions is always 
{ip\S\ip') instead of For instance, {ipml'ipn) 7^ ^mn between the eigenfunctions of ((1} 

because it is actually not physically meaningful, but {4'm\S\i'n) = {4'm\4'n) = ^mn is. (Please 
note that ip is used to denote the true wavefunction with nodal structure, and ip to denote 
pseudowavefunction, which are opposite in some papers.) 

H consists of the kinetic energy operator T, ionic local pseudopotential Vl, ionic nonlocal 
pseudopotential Vnl, Hartree potential Vh, and exchange-correlation potential Vxc, 

# = f + VL + V^L + Vk + Vxc. (2) 

The S operator is given by 

^ = i + E4l/5/)(A'l, (3) 
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where i = [Tim) is the angular momentum channel number, and / labels the ions. S contains 
contributions from all ions in the supercell, just as the total pseudopotential operator Vl + 
Vnl, which is the sum of pseudopotential operators of all ions. In above, the projector 
function /3/(x) = (x|/3/) of atom J's channel i is 

/?/(x)=A(x-X,), (4) 

where X/ is the ion position, and /3i(x) vanishes outside the pseudopotential cutoff. These 
projector functions appear in the nonlocal pseudopotential 

■1,3,1 

as well, where 

D'n = + / ^x(\/l(x) + \/h(x) + Vkc(x))Qj,(x). (6) 

The coefficients D^jf'^ are the unscreened scattering strengths, while the coefficients need 
to be self-consistently updated with the electron density 

p(x) = Y.\ + I /(^")' (7) 

in which is the Fermi-Dirac distribution. Qjj(x) is the charge augmentation function, 

i.e., the difference between the true wavefunction charge (interference) and the pseudocharge 
for selected channels, 

Q},(x) ^ (x)^/(x) - (x)^/(x), (8) 

which vanishes outside the cutoff. There is also 

ql^ ^ J cixgj,(x). (9) 

Terms in Eq. (|7j) are evaluated using two different grids, a sparse grid for the wavefunctions 
ipn and a dense grid for the augmentation functions Qjj(x). Ultrasoft pseudopotentials are 
thus fully specified by the functions Vl(x), /?/(x), Dji^\ and Qj^{x). Forces on ions and 
internal stress on the supercell can be derived analytically using linear response theory^^^. 

To extend the above ground-state USPP formalism to the time-dependent case, we note 
that the S operator in depends on the ionic positions {X/} only and not on the electronic 
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charge density. In the case that the ions are not moving, the following dynamical equations 
are equivalent: 

H{t)Mt) = thdtiSMt)) = s{thdtMt)), (10) 

whereby we have replaced the En in ((U) by the ihdt operator, and H{t) is updated using the 
time-dependent p(x, t). However when the ions are moving, 

ihdtS ^ s{ihdt) (11) 

with difference proportional to the ionic velocities. To resolve this ambiguity, we note that 
5* can be split as 

S = {S'/^U){U^S^/^), (12) 
where f/ is a unitary operator, f/f/^ = /, and we can rewrite (Q) as 

{U^S-'/')HiS-'/'U){U^S'/')^lj^ = e^{U^S'/')i;^. (13) 

Referring to the PAW formulation^, we can select U such that f/^S"^/^ is the PAW transfor- 
mation operator 

f/t5^/2 = f = l + J](|^/)-|^/))(/3/|: V^„ = m, (14) 
that maps the pseudowavefunction to the true wavefunction. So we can rewrite as, 

where H is then the true all-electron Hamiltonian (with core-level electrons frozen). In the 
all-electron TDDFT procedure, the above is replaced by the ihdt operator. It is thus 
clear that a physically meaningful TD-USPP equation in the case of moving ions should be 

or 

{u^s-'/^W^ = ihdt{{u^s'/')^n)- (17) 

In the equivalent PAW notation, it is simply, 

{P'r'H^n = indt{f^^). (18) 
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Or, in pseudized form amenable to numerical calculations, 

H^r. = thf\dt{f^n)) = in{f^f{dM + f\dtf)^n). (19) 

Differentiating (IHj), there is, 

a,f = (^^T^i/J/I + d"'.') - • x„ (20) 

and so we can define and calculate 

P = -ihf\dtf) = ^ • X/ (21) 

i,i 

operator, similar to analytic force calculation^. The TD-USPP / TD-PAW equation there- 
fore can be rearranged as, 

{H + P)^n = inS{dM, (22) 

with P proportional to the ionic velocities. It is basically the same as traditional TDDFT 
equation, but taking into account the moving spatial "gauge" due to ion motion. As such it 
can be used to model electron-phonon coupling^!, cluster dynamics under strong laser field^, 
etc., as long as the pseudopotential cores are not overlapping, and the core- level electrons 
are not excited. 

At each timestep, one should update p(x, t) as 

p(x,t) = 5^1 |^„(x,t)p + 5^Qj,(x)(^„(t)|/5/)(/5/|^„(t)) (23) 

n \ i,j,I J 

Note that while ipni^jt = 0) may be an eigenstate if we start from the ground-state wave- 
functions, ipni^it > 0) generally is no longer so with the external field turned on. n is 
therefore merely used as a label based on the initial state rather than an eigenstate label at 
t > 0. fn on the other hand always maintains its initial value, /„(t) = /^(O), for a particular 
simulation run. 

One may define projection operator tj belonging to atom J: 

i 

tj spatially has finite support, and so is 

dtj dij dil+ii) ,^ - 



Therefore in ()2ip is, 



P^ = -thf^^^' 



Oil 



= (l+tl)(l + t,)p-(l+tl)p(l + t,), (26) 

where p is the electron momentum operator. Unfortunately P'^ and therefore P are not 
Hermitian operators. This means that the numerical algorithm for integrating ()22j) may be 
different from the special case of immobile ions: 

H{t)i,n = thS{dtlPn). (27) 

Even if the same time-stepping algorithm is used, the error estimates would be different. In 
section III we discuss algorithms for integrating ipTj) only, and postpone detailed discussion 
of integration algorithm and error estimates for coupled ion-electron dynamics under 
USPP to a later paper. 

III. TIME-STEPPING ALGORITHMS FOR THE CASE OF IMMOBILE IONS 

In this section we focus on the important limiting case of (P7j) . where the ions are immobile 
or can be approximated as immobile. We may rewrite (jTrj) formally as 

S-'/'Hit)S-'/\S'/'^n) = ihd,{S'l^^^). (28) 
And so the time evolution of (j?fjl can be formally expressed as 



exp (^-'-j'^dt'S-^'^H{t')S-^'^^ 



^i/Vn(0), (29) 



with T the time-ordering operator. Algebraic expansions of different order are then per- 
formed on the above propagator, leading to various numerical time-stepping algorithms. 

A. First-order Implicit Euler Integration Scheme 

To first-order accuracy in time there are two well-known propagation algorithms, namely, 
the explicit (forward) Euler 

ihS — = H^n[^,t) (30) 
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and implicit (backward) Euler 

schemes. Although the explicit scheme ()30p is less expensive computationally, our test 
runs indicate that it always diverges numerically. The reason is that ()27p has poles on 
the imaginary axis, which are marginally outside of the stability domain {Iie{zAt) < 0) 
of the explicit algorithm. Therefore only the implicit algorithm can be used, which upon 
rearrangement is, 

V^„(t + At) = 5^„(t). (32) 

In the above, we still have the choice of whether to use H(t) or H(t + At). Since this 
is a first-order algorithm, neither choice would influence the order of the local truncation 
error. Through numerical tests we found that the implicit time differentiation in pip already 
imparts sufficient stability that the H{t + At) operator is not needed. Therefore we will 
solve 

^„(t + At) = S'^„(t) (33) 

at each timestep. Direct inversion turns out to be computationally infeasible in large- 
scale planewave calculations. We solve iteratively using matrix-free linear equation 
solvers such as the conjugate gradient method. Starting from the wavefunction of a previous 
timestep, we find that typically it takes about three to five conjugate gradient steps to achieve 
sufficiently convergent update. 

One serious drawback of this algorithm is that norm conservation of the wavefunction 

{Mt + ^t)\S\ijn{t + At)) = {tljn{t)\S\Mt)) (34) 

is not satisfied exactly, even if there is perfect floating-point operation accuracy. So one has 
to renormalize the wavefunction after several timesteps. 



S + ^HAt 



S + ^-H{t)At 



B. First-order Crank-Nicolson Integration Scheme 



We find the following Crank-Nicolson expansion2^*^2i54 of propagator 



SHr.it + At) = ^^£4S^^^^S^^^n(t) (35) 

^ ' l + l^S-'2H{t)S--2At ^ ' 



10 



stable enough for practical use. The norm of the wavefunction is conserved explicitly in the 
absence of roundoff errors, because of the spectral identity 



1 - l,S-^HS--^At 



1. 



(36) 



Therefore (j34|) is satisfied in an ideal numerical computation, and in practice one does not 
have to renormalize the wavefunctions in thousands of timesteps. 
Writing out the ()35|) expansion explicitly, we have: 



!/>„(« + At) 



(37) 



Similar to we solve Eq. (jTFj) using the conjugate gradient linear equations solver. This 
algorithm is still first-order because we use H{t), not {H{t) + H{t + At))/2, in ()37|) . In the 
limiting case of time- invariant charge density, p(x, t) = p(x, 0) and Hit + At) = H{t), the 
algorithm has second-order accuracy. This may happen if there is no external perturbation 
and we are simply testing whether the algorithm is stable in maintaining the eigenstate 
phase oscillation: ipnit) = '^n(0)e~*'^*, or in the case of propagating a test electron, which 
carries an infinitesimal charge and would not perturb H{t). 



C. Second-order Crank-Nicolson Integration Scheme 

We note that replacing H{t) by {H(t) + H{t + At))/2 in (jHH) would enhance the local 
truncation error to second order, while still maintaining norm conservation. In practice we of 
course do not know Hit + At) exactly, which depends on p{t + At) and therefore ipnit + At). 
However a sufficiently accurate estimate of p{t + At) can be obtained by running (|H7jl first 
for one step, from which we can get: 

p'{t + At) = p{t + At) + 0{At^), H'(t + At) = H{t + At) + 0{At^). (38) 



After this "predictor" step, we can solve: 
i{H{t) + H'it + At))At 



S + 



4h 



^nit + At) 



s- 



i{H{t)+H'{t + At))At 



^n(t), (39) 



to get the more accurate, second-order estimate for ipnit + At), that also satisfies (jMjl . 
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IV. OPTICAL ABSORPTION SPECTRA 



Calculating the optical absorption spectra of molecules, clusters and solids is one of 
the most important applications of TDDF T'^^i^^i^^i^^i^^i^^i^^i^^i^^i^''' . Since many experimen- 
tal and standard TDLDA results are available for comparison, we compute the spectra 
for sodium dimer (Na2) and benzene molecule {CqRq) to validate our direct time-stepping 
USPP-TDDFT scheme. 

We adopt the method by Bertsch et alM^ whereby an impulse electric field E(t) = 
eKk6(t)/e is applied to the system at t = 0, where k is unit vector and e is a small quantity. 
The system, which is at its ground state at t = 0^, would undergo transformation 

M^, t = 0+) = e^^'^-"^„(x, t = 0-), (40) 

for all its occupied electronic states, n = 1..N, at t = O"*". Note that the true, unpseudized 
wavefunctions should be used in ()4()|1 if theoretical rigor is to be maintained. 

One may then evolve {?/'„(x, t),n = l..A^} using a time stepper, with the total charge 
density p(x, t) updated at every step. The electric dipole moment d(t) is calculated as 

dW ^ e/.»xp(x..)x. (41) 

In a supercell calculation one needs to be careful to have a large enough vacuum region 
surrounding the molecule at the center, so no significant charge density can "spill over" the 
PBC boundary, thus causing a spurious discontinuity in d(t). 
The dipole strength tensor S(a;) can be computed by 

S(a;)k = m(u;) = lim - sinMe"^*' [d(t) - d(0)], (42) 

en-TT e,7^o e Jq 

where 7 is a small damping factor and rrie is the electron mass. In reality, the time integration 
is truncated at t{, and 7 should be chosen such that e~'''*f <^ 1. The merit of this and 
similar time-stepping approaches^U is that the entire spectrum can be obtained from just 
one calculation. 

For a molecule with no symmetry, one needs to carry out Eq. pn)) with subse- 
quent time integration for three independent k's: ki,k2,k3, and obtain three different 
mi(co'), m2{uj), 1113 (cj) on the right-hand side of Eq. (j42j) . One then solves the matrix equa- 
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tion: 



S(c^)[ki k2 ka] = [mi (w) 1112(0;) m3(cj)] S{lj) = [mi (u;) 1112(0;) 1113(0;)] [ki k2 ka] \ 

(43) 

S(o;) satisfies the Thomas-Reiche-Kuhn /-sum rule, 

POO 

N5ij = / diuSijiuj). (44) 
Jo 

For gas-phase systems where the orientation of the molecule or cluster is random, the 
isotropic average of S(o;) 

Sioo) ^ lTrS(o;) (45) 

may be calculated and plotted. 

In actual calculations employing norm-conserving pseudopotentials?-, the pseudo- 

wavefunctions '?/'„(x, t) are used in ()40|) instead of the true wavefunctions. And so the 

oscillator strength S(o;) obtained is not formally exact. However, the /-sum rule Eq. (I44|l 
is still satisfied exactly. With the USPP/PAW ioimalismM'^^'^^^, formally we should solve 

f V^„(x, t = 0+) = e^^'^-^f V^„(x, t = 0"), (46) 

using linear equation solver to get ?/'„(x, t = 0"*"), and then propagate '?/'„(x, t). However, 
for the present paper we skip this step, and replace ipn by ipn in (jinjl directly. This "quick- 
and-dirty fix" makes the oscillator strength not exact and also breaks the sum rule slightly. 
However, the peak positions are still correct. 

For the Na2 cluster, we actually use norm-conserving TM pseudopotential^'' for the Na 
atom, which is a special limiting case of our USPP-TDDFT code. The supercell is a tetrag- 

o 3 

onal box of 12 x 10 x 10 A and the Na2 cluster is along the x-direction with a bond length 
of 3.0 A. The planewave basis has a kinetic energy cutoff of 300 eV. The time integration 
is carried out for 10,000 steps with a timestep of At = 1.97 attoseconds, and e = 0.01/A, 
7 = 0.02eV^/^^. In the dipole strength plot (Fig. Q), the three peaks agree very well with 
TDLDA result from Octopus-^, and differ by ~ 0.4 eV from the experimental peaks^i^. In 
this case, the /-sum rule is verified to be satisfied to within 0.1% numerically. 

For the benzene molecule, ultrasoft pseudopotentials are used for both carbon and hy- 

o 3 

drogen atoms. The calculation is performed in a tetragonal box of 12.94 x 10 x 7 A with 
the benzene molecule placed on the x — y plane. The C-C bond length is 1.39 A and the 
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FIG. 1: Optical absorption spectra of Na2 cluster obtained from direct time-stepping TDLDA 
calculation using norm-conserving TM pseudopotential. The results should be compared with Fig. 
1 of Marques et ali^. 

C-H bond length is 1.1 A. The kinetic energy cutoff is 250 eV, e = 0.01/A, 7 = 0.1eVV^^ 
and the time integration is carried out for 5000 steps with a timestep of At = 2.37 at- 
toseconds. In the dipole strength function plot (Fig. 12), the peak at 6.95 eV represents 
the TT —>■ n* transition and the broad peak above 9 eV corresponds to the a —>■ a* transi- 
tion. The dipole strength function agrees very well with other TDLDA calculationa^SiS^ and 
experimenlj^. The slight difference is mostly due to our ad hoc approximation that ip^s 
instead of tpnS are used in (j40|) . The more formally rigorous implementation of the electric 
impulse perturbation, Eq. (j46|) . will be performed in future work. 

In this section we have verified the soundness of our time stepper with planewave basis 
through two examples of explicit electronic dynamics, where the charge density and effec- 
tive potential are updated at every timestep, employing both norm-conserving and ultrasoft 
pseudopotentials. This validation is important for the following non-perturbative propaga- 
tion of electrons in more complex systems. 
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FIG. 2: Optical absorption spectrum of benzene (CgHg) molecule. The results should be compared 
with Fig. 2 of Marques et al^ 

V. FERMI ELECTRON TRANSMISSION 

We first briefly review the setup of the Landauer transmission equation,i^>2^^ before 
performing an explicit TDDFT simulation. In its simplest form, two identical metallic leads 
(see Fig. Q) are connected to a device. The metallic lead is so narrow in y and z that only 
one channel (lowest quantum number in the y, z quantum well) needs to be considered. In 
the language of band structure, this means that one and only one branch of the ID band 
structure crosses the Fermi level Ep for kr^ > 0. Analogous to the universal density of states 
expression dN = 2Qdkxdkydkz/ {2ttY for 3D bulk metals, where Q is the volume and the 
factor of 2 accounts for up- and down-spins, the density of state of such ID system is simply 

d/V = (47) 
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In other words, the number of electrons per unit length with wave vector G {kx, + dkx) is 
just dkx/ 71- These electrons move with group velocity^: 

dE{kx) 



Vg 



fldkr 



(48) 



so there are [dkx / T^){dE{kx) / {fidkx)) = 2dE/h such electrons hitting the device from either 
side per unit time. 







left metal: Ep+edV/2 


device 


right metal: Ep-edV/2 



FIG. 3: Illustration of the Landauer transmission formalism. 

Under a small bias voltage dV, the Fermi level of the left lead is raised to Ep + edV/2, 
while that of the right lead drops to Ep — edV/2. The number of electrons hitting the device 
from the left with wave vector {kx, kx + dkx) is exactly equal to the number of electrons 
hitting the device from the right with wave vector {—kx,—kx — dkx), except in the small 
energy window {Ep — edV/2, Ep + edV/2), where the right has no electrons to balance against 
the left. Thus, a net number of 2{edV)/h electrons will attempt to cross from left and right, 
whose energies are very close to the original Ep. Some of them are scattered back by the 
device, and only a fraction of T G (0, 1] gets through. So the current they carry is: 



dl_ 

dV 



v=o 



'-^TiEp) 



(49) 



where 2e^/h = TTASlfiS = (12.906^)"!. 

Clearly, if the device is also of the same material and structure as the metallic leads, 
then T{Ep) should be 1, when we ignore electron-electron and electron-phonon scattering. 
This can be used as a sanity check of the code. For a nontrivial device however such as a 
molecular junction, T{Ep) would be smaller than 1, and would sensitively depend on the 
alignment of the molecular levels and Ep, as well as the overlap between these localized 
molecular states and the metallic states. 

Here we report two USPP-TDDFT case studies along the line of the above discussion. 
One is an infinite defect-free gold chain (Fig. I^a)). The other case uses gold chains as metal- 
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lie leads and eonneets them to a -S-C6H4-S- (benzene- (l,4)-ditliiolate, or BDT) molecular 
junction (Fig. Efb)). 




(b) 

FIG. 4: Atomistic configurations of our USPP-TDDFT simulations (Au: yellow, S: magenta, C: 
black, and H: white), (a) 12-atom Au chain. Bond length: Au-Au 2.88 A. (b) BDT (-S-C6H4-S-) 
junction connected to Au chain contacts. Bond lengths: Au-Au 2.88 A, Au-S 2.41 A, S-C 1.83 A, 
C-C 1.39 A, and C-H 1.1 A. 

In the semi-classical Landauer picture explained above, the metallic electrons are repre- 
sented by very wide Gaussian wavepacka^i moving along with the group velocity Vq, and 
with negligible rate of broadening compare to vq. Due to limitation of computational cost, 
we can only simulate rather small systems. In our experience with ID lithium and gold 
chains, a Gaussian envelop of 3-4 lattice constants in full width half maximum is suffi- 
cient to propagate at the Fermi velocity vcikp) with 100% transmissions and maintain its 
Gaussian-profile envelop with little broadening for several femto-seconds. 
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A. Fermi electron propagation in gold chain 



The ground-state electronic configurations of pure gold chains are calculated using the 
free USPP-DFT package DACAPO,ii'M with local density functional (LDA)M and 
planewave kinetic energy cutoff of 250 eV. The ultrasoft pseudopotential is generated using 
the free package uspp (ver. 7.3.3)^'^^^, with 5d, 6s, 6p, and auxiliary channels. Fig. Ufa) 
shows a chain of 12 Au atoms in a tetragonal supercell (34.56 x 12 x 12 A'^), with equal 
Au-Au bond length of 2.88 A. Theoretically, ID metal is always unstable against period- 
doubling Peierls distortion^i*^. However, the magnitude of the Peierls distortion is so small 
in the Au chain that room-temperature thermal fluctuations will readily erase its effect. For 
simplicity, we constrain the metallic chain to maintain single periodicity. Only the F-point 
wavefunctions are considered for the 12-atom configuration. 

The Fermi level Ep is found to be —6.65 eV, which is confirmed by a more accurate 
calculation of a one-Au-atom system with k-sampling (Fig. [S]). The Fermi state is doubly 
degenerate due to the time-inversion symmetry, corresponding to two Bloch wavefunctions 
of opposite wave vectors kp and —kp. 




0.1 0.2 0.3 m M 



FIG. 5: Band structure of a one- atom Au chain with 64 Monkliorst-PackSi k-sampling in the chain 
direction. The Fermi level, located at —6.65 eV, is marked as the dashed line. 
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From the F-point calculation, two energetically degenerate and real eigen-wavefunctions, 
ip+i'x) and ip-lx), are obtained. The complex traveling wavefunction is reconstructed as 

V^+(x) +2V^_(x) 



V2 



(50) 



The phase velocity of tpkp (x, t) computed from our TDLDA runs matches the Fermi frequency 
E-p/h. We use the integration scheme (jHTf) and a timestep of 2.37 attoseconds. 

We then calculate the Fermi electron group velocity wqI^f) by adding a perturbation 
modulation of 

V^fcp(x,t = 0) = V^fcp(x)(l + Asin(27rx/L)) (51) 

to the Fermi wavefunction ipk-p{^i where A is 0.02 and L is the x-length of the supercell. 
Fig. ini shows the electron density plot along two axes, x and t. From the line connecting the 
red-lobe edges, one can estimate the Fermi electron group velocity to be ~10.0 A/fs. The 
Fermi group velocity can also be obtained analytically from Eq. (j48j) at = k-p- A value 
of 10 A/fs is found according to Fig. El consistent with the TDLDA result. 



3.5 



0.5 



1 2 3 4 5 6 7 
Atom label 



9 10 11 12 



0.024 0.026 0.028 0.03 0.032 0.034 



FIG. 6: Evolution of modulated Fermi electron density in time along the chain direction. The 
electron density, in the unit of A ^, is an integral over the perpendicular y-z plane and normalized 
along the x direction, which is then color coded. 
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Lastly, the angular momentum projected densities of states are shown in Fig. [71 which 
indicate that the Fermi wavefunction mainly has s and px characteristics. 



T 1 1 1 1 1 1 I^Z^^ 




FIG. 7: Projected density of states of the 12-atom Au chain. 



B. Fermi electron transmission through Au-BDT-Au junction 

At small bias voltages, the electric conductance of a molecular junction (Fig. El^b)) is 
controlled by the transmission of Fermi electrons, as shown in Eq. (j49|) . In this section, we 
start from the Fermi electron wavefunction of a perfect ID gold chain (Fig. H^a)), and apply 
a Gaussian window centered at xq with a half width of cr, to obtain a localized wave pack 

V^,,(x,t = 0) = z^,,(x)G , (52) 

at the left lead. This localized Fermi electron wave pack is then propagated in real time by 
the TDLDA-USPP algorithm ()H7|1 with a timestep of 2.37 attoseconds, leaving from the left 
Au lead and traversing across the -S-C6H4-S- molecular junction (Fig. ^b)). While crossing 
the junction the electron will be scattered, after which we collect the electron density entering 
the right Au lead to compute the transmission probability T{E-p) literally. The calculation 
is performed in a tetragonal box (42.94 x 12 x 12 A'^) with a kinetic energy cutoff of 250 eV. 
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FIG. 8: Evolution of filtered wave package density in time along the chain direction. The electron 
density, in the unit of A ^, is a sum over the perpendicular y-z plane and normalized along the x 
direction. The normalized electron density is color coded by the absolute value. 

Fig. ISl shows the Fermi electron density evolution in x-t. A group velocity of 10 A/fs is 
obtained from the initial wave pack center trajectory, consistent with the perfect Au chain 
result. This free propagation lasts for about 0.8 fs, followed by a sharp density turnover that 
indicates the occurrence of strong electron scattering at the junction. A very small portion 
of the wave pack goes through the molecule. After about 1.7 fs, the reflected portion of the 
wave pack enters the right side of the supercell through PBC 

To separate the transmitted density from the reflected density as clearly as possible, we 
define and calculate the following cumulative charge on the right side 

PX' pLy 

R(x',t) = dx dy dzp{x^y, z,t), (53) 
Jxs Jo Jo 

where xs is the position of the right sulfur atom. R{x', t) is plotted in Fig. El for ten x'- 
positions starting from the right sulfur atom up to the right boundary L^. A shoulder can 
be seen in all 10 curves, at t = 1.5-2 fs, beyond which R{x',t) starts to rise sharply again, 
indicating that the reflected density has entered from the right boundary. Two solid curves 
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are highlighted in Fig. El The lower curve is at x' = xs + 7.2 A, which shows a clear 
transmission plateau of about 5%. The upper curve, which is for x' exactly at the right 
PBC boundary, shows R{x',t) ~ 7% at the shoulder. From these two curves, we estimate a 
transmission probability T{Ep) of 5-7%, which corresponds to a conductance of 4.0-5.6 fiS 
according to Eq. pHj) . This result from planewave TDLDA-USPP calculation is comparable 
to the transmission probability estimate of 10% from complex band structure calculations^*^ 
for one benzene linker (-C6H4-) without the sulfur atoms, and the non-equilibrium Green's 
function estimate of 5 /xSii for the similar system. 




0.5 1 1.5 2 2.5 3 3.5 

Time (fs) 



FIG. 9: R{x',t) versus time plot. Curves are measured in 10 different regions with different x' 
positions, which equally divide the region from the right S atom to the boundary on the right hand 
side. 



VI. SUMMARY 

In this work, we develop TDDFT based on Vanderbilt ultrasoft pseudopotentials and 
benchmark this USPP-TDDFT scheme by calculating optical absorption spectra, which 
agree with both experiments and other TDDFT calculations. We also demonstrate a new 
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approach to compute the electron conductance through single-molecule junction via wave 
pack propagation using TDDFT. The small conductance of 4.0-5.6 fiS is a result of our fixed 
band approximation, assuming the electron added was a small testing electron and therefore 
generated little disturbing effects of the incoming electrons on the electronic structure of the 
junction. This result is of the same order of magnitude as the results given by the Green's 
function and the complex band approaches, both requiring similar assumptions. 
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